FastEnz_comparison

FastEnz comparison

Setup shop

Set working environment

Input files you will need: .vcf file (data processing steps output, likely from ipyrad) and a population codes csv (column of pop codes in order of sequences in vcf - more info below).

setwd("/Users/mayroberts/Documents/Projects/FastEnztest/")

Load libraries

library(adegenet)

## Loading required package: ade4

## 
##    /// adegenet 2.1.11 is loaded ////////////
## 
##    > overview: '?adegenet'
##    > tutorials/doc/questions: 'adegenetWeb()' 
##    > bug reports/feature requests: adegenetIssues()

library(ape)
library(StAMPP)

## Loading required package: pegas

## Registered S3 method overwritten by 'pegas':
##   method      from
##   print.amova ade4

## 
## Attaching package: 'pegas'

## The following object is masked from 'package:ape':
## 
##     mst

## The following object is masked from 'package:ade4':
## 
##     amova

library(poppr)

## This is poppr version 2.9.8. To get started, type package?poppr
## OMP parallel support: unavailable

library(vcfR)

## 
##    *****       ***   vcfR   ***       *****
##    This is vcfR 1.15.0 
##      browseVignettes('vcfR') # Documentation
##      citation('vcfR') # Citation
##    *****       *****      *****       *****

## 
## Attaching package: 'vcfR'

## The following objects are masked from 'package:pegas':
## 
##     getINFO, write.vcf

library(devtools)

## Loading required package: usethis

library("BiocManager")

## Bioconductor version '3.21' is out-of-date; the current release version '3.22'
##   is available with R version '4.5'; see https://bioconductor.org/install

## 
## Attaching package: 'BiocManager'

## The following object is masked from 'package:devtools':
## 
##     install

library(dartRverse)

## **********************************************
## **** Welcome to dartRverse [Version 1.0.6] ****
## **********************************************

## ── Core dartRverse packages ────────────────────────────────────── dartRverse ──
## ✔ dartR.base 1.0.5     ✔ dartR.data 1.0.8
## ── Installed dartRverse packages   ─────────────────────────────── dartRverse ──
## ✔ dartR.captive 1.0.2     ✔ dartR.sim     0.70 
## ✔ dartR.popgen  1.0.0     ✔ dartR.spatial 1.0.3
## ── Not [yet] installed dartRverse packages ─────────────────────── dartRverse ──
## ✖ dartR.sexlinked

library(dartR.base)
library(dartR.captive)
library(dartR.data)
library(dartR.popgen)
library(dartR.sim)
library(dartR.spatial)
library(ggplot2)
library(dartR)

## Registered S3 methods overwritten by 'genetics':
##   method                      from    
##   [.haplotype                 pegas   
##   print.LD                    SNPassoc
##   print.LD.data.frame         SNPassoc
##   print.summary.LD.data.frame SNPassoc
##   summary.LD.data.frame       SNPassoc
## Registered S3 methods overwritten by 'dartR':
##   method      from      
##   cbind.dartR dartR.base
##   rbind.dartR dartR.base
## **** Welcome to the legacy dartR [Version 2.9.9.5 ] ****
## 
## ##############################
##  Be aware this the legacy version of dartR, which is no longer developed. Please consider using the new version of dartR, called 'dartRverse'. As before you can find the latest version on CRAN. dartRverse has the same functionality as dartR, but is more up-to-date and maintained. To reduce dependencies, and improve the overall performance, dartRverse consists of several packacages, that are maintained seperately to reduce dependencies . To install dartRverse please download the package dartRverse and follow the instructions therein. For more information please visit: https://github.com/green-striped-gecko/dartRverse.
## ##############################
## 
## Be aware that owing to CRAN requirements and compatibility reasons not all functions of the package may run after the basic installation, as some packages could still be missing. Hence for a most enjoyable experience we recommend to run the function 
## gl.install.vanilla.dartR()
## This installs all missing and required packages for your version of dartR. In case something fails during installation please refer to this tutorial: https://github.com/green-striped-gecko/dartR/wiki/Installation-tutorial.
## 
## For information how to cite dartR, please use:
## citation('dartR')
## Global verbosity is set to: 2
## 
## 
## **** Have fun using dartR! ****
## 
## Attaching package: 'dartR'
## 
## The following objects are masked from 'package:dartR.captive':
## 
##     gl.assign.grm, gl.assign.mahalanobis, gl.assign.pa, gl.assign.pca,
##     gl.filter.parent.offspring, gl.grm, gl.grm.network,
##     gl.plot.network, gl.report.parent.offspring, utils.assignment,
##     utils.assignment_2, utils.assignment_3, utils.assignment_4
## 
## The following objects are masked from 'package:dartR.spatial':
## 
##     gl.costdistances, gl.genleastcost, gl.ibd, gl.spatial.autoCorr,
##     gl2shp, utils.spautocor
## 
## The following objects are masked from 'package:dartR.popgen':
## 
##     gl.blast, gl.collapse, gl.evanno, gl.ld.distance, gl.ld.haplotype,
##     gl.LDNe, gl.map.structure, gl.nhybrids, gl.outflank,
##     gl.plot.faststructure, gl.plot.structure, gl.run.faststructure,
##     gl.run.structure, gl.sfs, utils.outflank,
##     utils.outflank.MakeDiploidFSTMat, utils.outflank.plotter,
##     utils.structure.evanno, utils.structure.genind2gtypes,
##     utils.structure.run
## 
## The following objects are masked from 'package:dartR.sim':
## 
##     gl.diagnostics.sim, gl.sim.create_dispersal, gl.sim.emigration,
##     gl.sim.ind, gl.sim.mutate, gl.sim.offspring, gl.sim.WF.run,
##     gl.sim.WF.table
## 
## The following objects are masked from 'package:dartR.base':
## 
##     gi2gl, gl.alf, gl.allele.freq, gl.amova, gl.check.verbosity,
##     gl.check.wd, gl.colors, gl.compliance.check, gl.define.pop,
##     gl.diagnostics.hwe, gl.dist.ind, gl.dist.pop, gl.drop.ind,
##     gl.drop.loc, gl.drop.pop, gl.edit.recode.ind, gl.edit.recode.pop,
##     gl.fdsim, gl.filter.allna, gl.filter.callrate, gl.filter.hamming,
##     gl.filter.heterozygosity, gl.filter.hwe, gl.filter.locmetric,
##     gl.filter.maf, gl.filter.monomorphs, gl.filter.overshoot,
##     gl.filter.pa, gl.filter.rdepth, gl.filter.reproducibility,
##     gl.filter.secondaries, gl.filter.taglength, gl.fixed.diff,
##     gl.fst.pop, gl.He, gl.Ho, gl.hwe.pop, gl.impute, gl.join,
##     gl.keep.ind, gl.keep.loc, gl.keep.pop, gl.load, gl.make.recode.ind,
##     gl.make.recode.pop, gl.map.interactive, gl.merge.pop, gl.pcoa,
##     gl.pcoa.plot, gl.plot.heatmap, gl.print.history, gl.propShared,
##     gl.read.csv, gl.read.dart, gl.read.fasta, gl.read.silicodart,
##     gl.read.vcf, gl.reassign.pop, gl.recalc.metrics, gl.recode.ind,
##     gl.recode.pop, gl.rename.pop, gl.report.bases, gl.report.callrate,
##     gl.report.diversity, gl.report.fstat, gl.report.hamming,
##     gl.report.heterozygosity, gl.report.hwe, gl.report.ld.map,
##     gl.report.locmetric, gl.report.maf, gl.report.monomorphs,
##     gl.report.overshoot, gl.report.pa, gl.report.rdepth,
##     gl.report.replicates, gl.report.reproducibility,
##     gl.report.secondaries, gl.report.taglength, gl.sample, gl.save,
##     gl.select.colors, gl.select.shapes, gl.set.verbosity, gl.smearplot,
##     gl.sort, gl.subsample.loci, gl.test.heterozygosity, gl.tree.nj,
##     gl.write.csv, gl2bayescan, gl2bpp, gl2demerelate, gl2eigenstrat,
##     gl2fasta, gl2faststructure, gl2gds, gl2genalex, gl2genepop,
##     gl2geno, gl2gi, gl2hiphop, gl2phylip, gl2plink, gl2related, gl2sa,
##     gl2structure, gl2treemix, gl2vcf, utils.basic.stats,
##     utils.check.datatype, utils.jackknife, utils.plot.save
## 
## The following objects are masked from 'package:dartR.data':
## 
##     bandicoot.gl, possums.gl, testset.gl, testset.gs

Import data:

Here we import the VCF file from ipyrad and then convert to genlight object. Read in the VCF file:

FastEnzTest.vcf <- read.vcfR("/Users/mayroberts/Documents/Projects/FastEnzTest/ipyrad/FastEnz_test_outfiles/FastEnz_test.vcf")

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 10
##   header_line: 11
##   variant count: 5578
##   column count: 27
## Meta line 10 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 5578
##   Character matrix gt cols: 27
##   skip: 0
##   nrows: 5578
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant 4000Processed variant 5000Processed variant: 5578
## All variants processed

FastEnzvsSNI.vcf <- read.vcfR("/Users/mayroberts/Documents/Projects/FastEnztest/ipyrad/FastEnz_only_outfiles/FastEnz_only.vcf")

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 10
##   header_line: 11
##   variant count: 273
##   column count: 21
## Meta line 10 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 273
##   Character matrix gt cols: 21
##   skip: 0
##   nrows: 273
##   row_num: 0
## Processed variant: 273
## All variants processed

OrigEnzvsSNI.vcf <- read.vcfR("/Users/mayroberts/Documents/Projects/FastEnztest/ipyrad/OrigEnz_outfiles/OrigEnz.vcf")

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 10
##   header_line: 11
##   variant count: 209
##   column count: 21
## Meta line 10 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 209
##   Character matrix gt cols: 21
##   skip: 0
##   nrows: 209
##   row_num: 0
## Processed variant: 209
## All variants processed

FastEnzOrig_only.vcf <- read.vcfR("/Users/mayroberts/Documents/Projects/FastEnztest/ipyrad/FastEnzOrig_only_outfiles/FastEnzOrig_only.vcf")

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 10
##   header_line: 11
##   variant count: 32602
##   column count: 21
## Meta line 10 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 32602
##   Character matrix gt cols: 21
##   skip: 0
##   nrows: 32602
##   row_num: 0
## Processed variant 1000Processed variant 2000Processed variant 3000Processed variant 4000Processed variant 5000Processed variant 6000Processed variant 7000Processed variant 8000Processed variant 9000Processed variant 10000Processed variant 11000Processed variant 12000Processed variant 13000Processed variant 14000Processed variant 15000Processed variant 16000Processed variant 17000Processed variant 18000Processed variant 19000Processed variant 20000Processed variant 21000Processed variant 22000Processed variant 23000Processed variant 24000Processed variant 25000Processed variant 26000Processed variant 27000Processed variant 28000Processed variant 29000Processed variant 30000Processed variant 31000Processed variant 32000Processed variant: 32602
## All variants processed

FastEnzTest.vcf

## ***** Object of Class vcfR *****
## 18 samples
## 2232 CHROMs
## 5,578 variants
## Object size: 2.4 Mb
## 0 percent missing data
## *****        *****         *****

OrigEnzvsSNI.vcf

## ***** Object of Class vcfR *****
## 12 samples
## 107 CHROMs
## 209 variants
## Object size: 0.1 Mb
## 0 percent missing data
## *****        *****         *****

FastEnzvsSNI.vcf

## ***** Object of Class vcfR *****
## 12 samples
## 125 CHROMs
## 273 variants
## Object size: 0.2 Mb
## 0 percent missing data
## *****        *****         *****

FastEnzOrig_only.vcf

## ***** Object of Class vcfR *****
## 12 samples
## 15280 CHROMs
## 32,602 variants
## Object size: 9.5 Mb
## 0 percent missing data
## *****        *****         *****

Convert the VCF file to genlight object:
Name it something easy to find as your genlight object, you will use it many times in the downstream analyses.

gl_FastEnzTest <- vcfR2genlight(FastEnzTest.vcf)

## Warning in vcfR2genlight(FastEnzTest.vcf): Found 62 loci with more than two alleles.
## Objects of class genlight only support loci with two alleles.
## 62 loci will be omitted from the genlight object.

gl_Orig <- vcfR2genlight(OrigEnzvsSNI.vcf)

## Warning in vcfR2genlight(OrigEnzvsSNI.vcf): Found 1 loci with more than two alleles.
## Objects of class genlight only support loci with two alleles.
## 1 loci will be omitted from the genlight object.

gl_Fast <- vcfR2genlight(FastEnzvsSNI.vcf)

## Warning in vcfR2genlight(FastEnzvsSNI.vcf): Found 1 loci with more than two alleles.
## Objects of class genlight only support loci with two alleles.
## 1 loci will be omitted from the genlight object.

gl_both <- vcfR2genlight(FastEnzOrig_only.vcf)

## Warning in vcfR2genlight(FastEnzOrig_only.vcf): Found 215 loci with more than two alleles.
## Objects of class genlight only support loci with two alleles.
## 215 loci will be omitted from the genlight object.

Import population codes

You will need to attach population fields to do any analyses by population. This is done by reading in a CSV file that is a single column of population codes/names and adding it as a new column to your genlight object.

Making your popcodes.csv The column of population codes should be in the same order as the order samples in your output from IPYRAD. Get this from the _stats.txt file in the _outfiles in your ipyrad directory and make the list in excel and a text editor and save as a .csv. It’s important that the order of the pop_codes matches the order of the sample names in the vcf file. Also, no spaces in site names.

Read in population file:

#this file has pop codes in column1 (V1)
pop.data <- read.table("/Users/mayroberts/Documents/Projects/FastEnztest/pop_file.csv")
pop.dataF <- read.table("/Users/mayroberts/Documents/Projects/FastEnztest/pop_file_F.csv")
pop.dataO <- read.table("/Users/mayroberts/Documents/Projects/FastEnztest/pop_file_O.csv")
pop.dataB <- read.table("/Users/mayroberts/Documents/Projects/FastEnztest/pop_file_both.csv")
str(pop.data)

## 'data.frame':    18 obs. of  2 variables:
##  $ V1: chr  "Lomatium_insulare_10727_GUZA" "Lomatium_insulare_10728_GUZA" "Lomatium_insulare_10729_GUZA" "Lomatium_insulare_10730_GUPS" ...
##  $ V2: chr  "GUorig" "GUorig" "GUorig" "GUorig" ...

str(pop.dataF)

## 'data.frame':    12 obs. of  2 variables:
##  $ V1: chr  "Lomatium_insulare_9900_SNI1" "Lomatium_insulare_9901_SNI1" "Lomatium_insulare_9902_SNI1" "Lomatium_insulare_9917_SNI3" ...
##  $ V2: chr  "SNI" "SNI" "SNI" "SNI" ...

str(pop.dataO)

## 'data.frame':    12 obs. of  2 variables:
##  $ V1: chr  "Lomatium_insulare_10727_GUZA" "Lomatium_insulare_10728_GUZA" "Lomatium_insulare_10729_GUZA" "Lomatium_insulare_10730_GUPS" ...
##  $ V2: chr  "GUorig" "GUorig" "GUorig" "GUorig" ...

str(pop.dataB)

## 'data.frame':    12 obs. of  1 variable:
##  $ V1: chr  "Lomatium_insulare_10727_GUZA\tGUorig" "Lomatium_insulare_10728_GUZA\tGUorig" "Lomatium_insulare_10729_GUZA\tGUorig" "Lomatium_insulare_10730_GUPS\tGUorig" ...

Add population codes to your genlight object

pop(gl_FastEnzTest) <- pop.data$V2
pop(gl_Fast) <- pop.dataF$V2
pop(gl_Orig) <- pop.dataO$V2
pop(gl_both) <- pop.dataB$V2
head(gl_FastEnzTest)

##  /// GENLIGHT OBJECT /////////
## 
##  // 6 genotypes,  5,516 binary SNPs, size: 643.1 Kb
##  10208 (30.84 %) missing data
## 
##  // Basic content
##    @gen: list of 6 SNPbin
## 
##  // Optional content
##    @ind.names:  6 individual labels
##    @loc.names:  5516 locus labels
##    @chromosome: factor storing chromosomes of the SNPs
##    @position: integer storing positions of the SNPs
##    @pop: population of each individual (group size range: 6-6)
##    @other: a list containing: elements without names

head(gl_Fast)

##  /// GENLIGHT OBJECT /////////
## 
##  // 6 genotypes,  272 binary SNPs, size: 42.4 Kb
##  88 (5.39 %) missing data
## 
##  // Basic content
##    @gen: list of 6 SNPbin
## 
##  // Optional content
##    @ind.names:  6 individual labels
##    @loc.names:  272 locus labels
##    @chromosome: factor storing chromosomes of the SNPs
##    @position: integer storing positions of the SNPs
##    @pop: population of each individual (group size range: 6-6)
##    @other: a list containing: elements without names

head(gl_Orig)

##  /// GENLIGHT OBJECT /////////
## 
##  // 6 genotypes,  208 binary SNPs, size: 36.2 Kb
##  100 (8.01 %) missing data
## 
##  // Basic content
##    @gen: list of 6 SNPbin
## 
##  // Optional content
##    @ind.names:  6 individual labels
##    @loc.names:  208 locus labels
##    @chromosome: factor storing chromosomes of the SNPs
##    @position: integer storing positions of the SNPs
##    @pop: population of each individual (group size range: 6-6)
##    @other: a list containing: elements without names

head(gl_both)

##  /// GENLIGHT OBJECT /////////
## 
##  // 6 genotypes,  32,387 binary SNPs, size: 4 Mb
##  115321 (59.35 %) missing data
## 
##  // Basic content
##    @gen: list of 6 SNPbin
## 
##  // Optional content
##    @ind.names:  6 individual labels
##    @loc.names:  32387 locus labels
##    @chromosome: factor storing chromosomes of the SNPs
##    @position: integer storing positions of the SNPs
##    @other: a list containing: elements without names

Quick check of data:

Genlight compliance check

This is run by dartR to make sure that the genlight object conforms to dartR requirements, if it doesn’t, it will rectify it. Can take a minute.

gl_FastEnzTest <- gl.compliance.check(gl_FastEnzTest)

## Starting gl.compliance.check 
##   Processing genlight object with SNP data
##   The slot loc.all, which stores allele name for each locus, is empty. 
## Creating a dummy variable (A/C) to insert in this slot. 
##   Checking coding of SNPs
##     SNP data scored NA, 0, 1 or 2 confirmed
##   Checking locus metrics and flags
##   Recalculating locus metrics
##   Checking for monomorphic loci
##     No monomorphic loci detected
##   Checking for loci with all missing data
##     No loci with all missing data detected
##   Checking whether individual names are unique.
##   Checking for individual metrics
##   Warning: Creating a slot for individual metrics
##   Checking for population assignments
##     Population assignments confirmed
##   Spelling of coordinates checked and changed if necessary to 
##             lat/lon
## Completed: gl.compliance.check

gl_Fast <- gl.compliance.check(gl_Fast)

## Starting gl.compliance.check 
##   Processing genlight object with SNP data
##   The slot loc.all, which stores allele name for each locus, is empty. 
## Creating a dummy variable (A/C) to insert in this slot. 
##   Checking coding of SNPs
##     SNP data scored NA, 0, 1 or 2 confirmed
##   Checking locus metrics and flags
##   Recalculating locus metrics
##   Checking for monomorphic loci
##     No monomorphic loci detected
##   Checking for loci with all missing data
##     No loci with all missing data detected
##   Checking whether individual names are unique.
##   Checking for individual metrics
##   Warning: Creating a slot for individual metrics
##   Checking for population assignments
##     Population assignments confirmed
##   Spelling of coordinates checked and changed if necessary to 
##             lat/lon
## Completed: gl.compliance.check

gl_Orig <- gl.compliance.check(gl_Orig)

## Starting gl.compliance.check 
##   Processing genlight object with SNP data
##   The slot loc.all, which stores allele name for each locus, is empty. 
## Creating a dummy variable (A/C) to insert in this slot. 
##   Checking coding of SNPs
##     SNP data scored NA, 0, 1 or 2 confirmed
##   Checking locus metrics and flags
##   Recalculating locus metrics
##   Checking for monomorphic loci
##     No monomorphic loci detected
##   Checking for loci with all missing data
##     No loci with all missing data detected
##   Checking whether individual names are unique.
##   Checking for individual metrics
##   Warning: Creating a slot for individual metrics
##   Checking for population assignments
##     Population assignments confirmed
##   Spelling of coordinates checked and changed if necessary to 
##             lat/lon
## Completed: gl.compliance.check

gl_both <- gl.compliance.check(gl_both)

## Starting gl.compliance.check 
##   Processing genlight object with SNP data
##   The slot loc.all, which stores allele name for each locus, is empty. 
## Creating a dummy variable (A/C) to insert in this slot. 
##   Checking coding of SNPs
##     SNP data scored NA, 0, 1 or 2 confirmed
##   Checking locus metrics and flags
##   Recalculating locus metrics
##   Checking for monomorphic loci
##     No monomorphic loci detected
##   Checking for loci with all missing data
##     No loci with all missing data detected
##   Checking whether individual names are unique.
##   Checking for individual metrics
##   Warning: Creating a slot for individual metrics
##   Checking for population assignments
##   Population assignments not detected, individuals assigned
##                     to a single population labelled 'pop1'
##   Spelling of coordinates checked and changed if necessary to 
##             lat/lon
## Completed: gl.compliance.check

N per island

This counts samples by population and makes and object calls it popvsN

popvsn <-table(pop(gl_FastEnzTest))
popvsn # table of populations and number of samples in them 

## 
## GUFast GUorig    SNI 
##      6      6      6

barplot(popvsn, las=2, cex.names = .6) #Barplot of number of samples per population

popvsn <-table(pop(gl_Fast))
popvsn # table of populations and number of samples in them 

## 
## GUFast    SNI 
##      6      6

barplot(popvsn, las=2, cex.names = .6) #Barplot of number of samples per population

popvsn <-table(pop(gl_Orig))
popvsn # table of populations and number of samples in them 

## 
## GUorig    SNI 
##      6      6

barplot(popvsn, las=2, cex.names = .6) #Barplot of number of samples per population

### Assess missing data
This produces a plot where each sample is in a row down the y axis and loci are across the x. Can take a minute or two depending on number of samples and seqs.

glPlot(gl_FastEnzTest, posi="topleft")
title("FastEnzTest_all glPlot")

glPlot(gl_both, posi="topleft")
title("FastEnz and OrigEnz glPlot")

# Convert to matrix
mat <- as.matrix(gl_FastEnzTest)

#Adjust margins: bottom, left, top, right
# Increase left margin (second number) for long sample names
par(mar = c(5, 12, 4, 2))  

# Plot the heatmap
image(t(mat), axes = FALSE, col = c("navy", "gold", "red"))

# Add y-axis labels
axis(2, at = seq(0, 1, length.out = nInd(gl_FastEnzTest)),
     labels = indNames(gl_FastEnzTest), las = 2, cex.axis = 0.5)

# Add box around the plot
box()

# Add legend
legend("topright", legend = c("0 = homo ref", "1 = het", "2 = homo alt"),
       fill = c("navy", "gold", "red"), cex = 0.7, bty = "n")

Ordered so that pairs of samples are next to each other

########
# Custom order of sample names
my_order <- c("Lomatium_insulare_10727_GUZA", "Lomatium_insulare_fastEnz_10727_GUZA", "Lomatium_insulare_10728_GUZA", "Lomatium_insulare_fastEnz_10728_GUZA", "Lomatium_insulare_10729_GUZA", "Lomatium_insulare_fastEnz_10729_GUZA", "Lomatium_insulare_10730_GUPS", "Lomatium_insulare_fastEnz_10730_GUPS", "Lomatium_insulare_10731_GUPS", "Lomatium_insulare_fastEnz_10731_GUPS", "Lomatium_insulare_10737_GUPS", "Lomatium_insulare_fastEnz_10737_GUPS","Lomatium_insulare_9900_SNI1", "Lomatium_insulare_9901_SNI1","Lomatium_insulare_9902_SNI1", "Lomatium_insulare_9917_SNI3", "Lomatium_insulare_9918_SNI3", "Lomatium_insulare_9919_SNI3")
# Reorder the genlight object
gl_FastEnzTest_ordered <- gl_FastEnzTest[match(my_order, indNames(gl_FastEnzTest))]

# Convert to matrix
mat <- as.matrix(gl_FastEnzTest_ordered)

#Adjust margins: bottom, left, top, right
# Increase left margin (second number) for long sample names
par(mar = c(5, 12, 4, 2))  

# Plot the heatmap
image(t(mat), axes = FALSE, col = c("navy", "gold", "red"))

# Add y-axis labels
axis(2, at = seq(0, 1, length.out = nInd(gl_FastEnzTest_ordered)),
     labels = indNames(gl_FastEnzTest_ordered), las = 2, cex.axis = 0.5)

# Add box around the plot
box()

# Add legend
legend("topright", legend = c("0 = homo ref", "1 = het", "2 = homo alt"),
       fill = c("navy", "gold", "red"), cex = 0.7, bty = "n")

## Neighbor joining tree ### NJ tree - colors by site I recommend running this in the console directly to get saveable plots in the “plots” window. Make sure the plots window is expanded.

library(ape)
#Distance matrix as tre object
tre <- njs(dist(as.matrix(gl_FastEnzTest)))
treF <- njs(dist(as.matrix(gl_Fast)))
treO <- njs(dist(as.matrix(gl_Orig)))
treB <- njs(dist(as.matrix(gl_both)))
#

Plot FastEnz + Orig + SNI together

# Define color mapping
pop_levels <- unique(gl_FastEnzTest$pop)  # e.g. c("SNI","GUorig","GUFast")
palette3 <- c("#0072B2","#009E73","#E69F00")  # length must match number of populations
names(palette3) <- pop_levels  # name colors by population

# Map tip colors
tip_colors <- palette3[as.character(gl_FastEnzTest$pop)]

# Plot tree
plot(tre, type="unrooted", show.tip.lab=TRUE, cex=1, font=1, tip.col=tip_colors)

# Add legend (line colors)
legend("topleft", legend = pop_levels,
       col = palette3, cex=1, fill = palette3, title = "Groups by color")
title("FastEnzyme, original enzyme comparison")

# Define color mapping
pop_levels <- unique(gl_both$pop)  # e.g. c("SNI","GUorig","GUFast")
palette3 <- c("#0072B2","#009E73","#E69F00")  # length must match number of populations
names(palette3) <- pop_levels  # name colors by population

# Map tip colors
tip_colors <- palette3[as.character(gl_both$pop)]

# Plot tree
plot(treB, type="unrooted", show.tip.lab=TRUE, cex=1, font=1, tip.col=tip_colors)

# Add legend (line colors)
legend("topleft", legend = pop_levels,
       col = palette3, cex=1, fill = palette3, title = "Groups by color")
title("FastEnzyme and original enzyme comparison")

#pdf("Fast_Orig_trees_large.pdf", width = 14, height = 7)  # in inches

# Set up 1 row, 2 columns
par(mfrow = c(1, 2))

## ---- Plot 1: Fast + SNI ----
pop_levels <- unique(gl_Fast$pop)
palette3 <- c("#E69F00","#009E73")
names(palette3) <- pop_levels
tip_colors <- palette3[as.character(gl_Fast$pop)]

plot(treF, type = "unrooted", show.tip.lab = TRUE, cex = .6, font = 1, tip.col = tip_colors)
legend("topleft", legend = pop_levels, fill = palette3, cex = 1, title = "Groups by color")
title("FastEnzyme vs SNI comparison")

## ---- Plot 2: Orig + SNI ----
pop_levels <- unique(gl_Orig$pop)
palette3 <- c("#0072B2","#E69F00")
names(palette3) <- pop_levels
tip_colors <- palette3[as.character(gl_Orig$pop)]

plot(treO, type = "unrooted", show.tip.lab = TRUE, cex = .6, font = 1, tip.col = tip_colors)
legend("topleft", legend = pop_levels, fill = palette3, cex = 1, title = "Groups by color")
title("OriginalEnzymes vs SNI comparison")

# Reset layout afterward
par(mfrow = c(1, 1))

#dev.off()

PCA

ran run_plink.sh in terminal with LOMINS_GUAD.vcf file

then to visualize:

library(ggplot2)
library(RColorBrewer)

set input variables

read in the output from plink

eigenval=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/fastEnz.eigenval")
pca=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/fastEnz.eigenvec")

eigenval_F=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/fastEnz_F.eigenval")
pca_F=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/fastEnz_F.eigenvec")

eigenval_O=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/OrigEnz.eigenval")
pca_O=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/OrigEnz.eigenvec")

eigenval_B=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/FastEnzOrig_only.eigenval")
pca_B=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/FastEnzOrig_only.eigenvec")

Import lists of samples by how you want them colored in your plot. Single column of sample name with no .filetype at ending ### FastEnz+OriginalEnz+SNI PCA plot

## FastEnz+OriginalEnz+SNI PCA plot
FastEnz=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/FastEnz.txt")
OrigEnz=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/OrigEnz.txt")
SNI=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/SNI.txt")
pca$Sites <- NA  # create empty column first

pca$Sites=ifelse(pca$V1%in%FastEnz$V1,"FastEnz",pca$Sites)
pca$Sites=ifelse(pca$V1%in%OrigEnz$V1,"OrigEnz",pca$Sites)
pca$Sites=ifelse(pca$V1%in%SNI$V1,"SNI",pca$Sites)

cols <- c("#0072B2","#009E73","#E69F00")

### Global plotting controls
dataset_name <- "EnzTest"
jitter_width  <- 0.05
jitter_height <- 0.05
alpha_val     <- 0.6
point_size    <- 4

### Function to plot any two PCs and save to file
plot_pca <- function(df, eigenval, pc_x, pc_y, cols, outdir = ".", format = "pdf") {
  
  # Build the ggplot object
  p <- ggplot(df, aes_string(x = paste0("V", pc_x+2), y = paste0("V", pc_y+2), colour = "Sites")) +
    geom_jitter(width = jitter_width, height = jitter_height, alpha = alpha_val, size = point_size) +
    scale_color_manual(values = cols) +
    geom_hline(yintercept = 0) +
    geom_vline(xintercept = 0) +
    theme_bw() +
    labs(
      x = paste0("PC", pc_x, " (", round(100 * eigenval[pc_x,1] / sum(eigenval[,1]), 2), "%)"),
      y = paste0("PC", pc_y, " (", round(100 * eigenval[pc_y,1] / sum(eigenval[,1]), 2), "%)")
    ) +
    theme(text = element_text(size = 10))
  
  # Construct filename
  filename <- file.path(outdir, paste0(dataset_name,"_PC", pc_x, "_PC", pc_y, ".", format))
  
  # Save depending on chosen format
  if (format == "pdf") {
    ggsave(filename, plot = p, device = "pdf", width = 6, height = 5)
  } else if (format == "png") {
    ggsave(filename, plot = p, device = "png", width = 6, height = 5, dpi = 300)
  } else {
    warning("Unsupported format. Use 'pdf' or 'png'.")
  }
  
  # Also return plot to R session
  return(p)
}

### plots
# Save to PDF (default) in current working directory
plot_pca(pca, eigenval, 1, 2, cols)    # PC1 vs PC2

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot_pca(pca, eigenval, 1, 3, cols)    # PC1 vs PC3

#plot_pca(pca, eigenval, 1, 4, cols)    # PC1 vs PC4
#plot_pca(pca, eigenval, 1, 5, cols)    # PC1 vs PC5
#plot_pca(pca, eigenval, 1, 6, cols)    # PC1 vs PC6
plot_pca(pca, eigenval, 2, 3, cols)    # PC2 vs PC3

plot_pca(pca, eigenval, 2, 4, cols)    # PC2 vs PC4

#plot_pca(pca, eigenval, 2, 5, cols)    # PC2 vs PC5
#plot_pca(pca, eigenval, 2, 6, cols)    # PC2 vs PC6
plot_pca(pca, eigenval, 3, 4, cols)    # PC3 vs PC4

#plot_pca(pca, eigenval, 3, 5, cols)    # PC3 vs PC5
#plot_pca(pca, eigenval, 3, 6, cols)    # PC3 vs PC6


# Save to PNG in a specific folder
#plot_pca(pca, eigenval, 1, 2, cols, outdir = "/Users/mayroberts/Documents/Projects/LOMINS/Analysis/plink/PCA_locin25psampls/plots", format = "png")

FastEnz+OriginalEnz PCA plot

## FastEnz+OriginalEnz+SNI PCA plot
FastEnz=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/FastEnz.txt")
OrigEnz=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/OrigEnz.txt")
#SNI=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/SNI.txt")
pca_B$Sites <- NA  # create empty column first

pca_B$Sites=ifelse(pca_B$V1%in%FastEnz$V1,"FastEnz",pca_B$Sites)
pca_B$Sites=ifelse(pca_B$V1%in%OrigEnz$V1,"OrigEnz",pca_B$Sites)
pca_B$Sites=ifelse(pca_B$V1%in%SNI$V1,"SNI",pca_B$Sites)

cols <- c("#0072B2","#009E73","#E69F00")

### Global plotting controls
dataset_name <- "EnzTest"
jitter_width  <- 0.05
jitter_height <- 0.05
alpha_val     <- 0.6
point_size    <- 4

### Function to plot any two PCs and save to file
plot_pca <- function(df, eigenval, pc_x, pc_y, cols, outdir = ".", format = "pdf") {
  
  # Build the ggplot object
  p <- ggplot(df, aes_string(x = paste0("V", pc_x+2), y = paste0("V", pc_y+2), colour = "Sites")) +
    geom_jitter(width = jitter_width, height = jitter_height, alpha = alpha_val, size = point_size) +
    scale_color_manual(values = cols) +
    geom_hline(yintercept = 0) +
    geom_vline(xintercept = 0) +
    theme_bw() +
    labs(
      x = paste0("PC", pc_x, " (", round(100 * eigenval[pc_x,1] / sum(eigenval[,1]), 2), "%)"),
      y = paste0("PC", pc_y, " (", round(100 * eigenval[pc_y,1] / sum(eigenval[,1]), 2), "%)")
    ) +
    theme(text = element_text(size = 10))
  
  # Construct filename
  filename <- file.path(outdir, paste0(dataset_name,"_PC", pc_x, "_PC", pc_y, ".", format))
  
  # Save depending on chosen format
  if (format == "pdf") {
    ggsave(filename, plot = p, device = "pdf", width = 6, height = 5)
  } else if (format == "png") {
    ggsave(filename, plot = p, device = "png", width = 6, height = 5, dpi = 300)
  } else {
    warning("Unsupported format. Use 'pdf' or 'png'.")
  }
  
  # Also return plot to R session
  return(p)
}

### plots
# Save to PDF (default) in current working directory
plot_pca(pca_B, eigenval, 1, 2, cols)    # PC1 vs PC2

plot_pca(pca_B, eigenval, 1, 3, cols)    # PC1 vs PC3

#plot_pca(pca, eigenval, 1, 4, cols)    # PC1 vs PC4
#plot_pca(pca, eigenval, 1, 5, cols)    # PC1 vs PC5
#plot_pca(pca, eigenval, 1, 6, cols)    # PC1 vs PC6
plot_pca(pca_B, eigenval, 2, 3, cols)    # PC2 vs PC3

plot_pca(pca_B, eigenval, 2, 4, cols)    # PC2 vs PC4

#plot_pca(pca, eigenval, 2, 5, cols)    # PC2 vs PC5
#plot_pca(pca, eigenval, 2, 6, cols)    # PC2 vs PC6
plot_pca(pca_B, eigenval, 3, 4, cols)    # PC3 vs PC4

#plot_pca(pca, eigenval, 3, 5, cols)    # PC3 vs PC5
#plot_pca(pca, eigenval, 3, 6, cols)    # PC3 vs PC6


# Save to PNG in a specific folder
#plot_pca(pca, eigenval, 1, 2, cols, outdir = "/Users/mayroberts/Documents/Projects/LOMINS/Analysis/plink/PCA_locin25psampls/plots", format = "png")

FastEnz only + SNI PCA plot

## FastEnz only + SNI PCA plot
FastEnz=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/FastEnz.txt")
SNI=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/SNI.txt")
pca_F$Sites <- NA  # create empty column first

pca_F$Sites=ifelse(pca_F$V1%in%FastEnz$V1,"FastEnz",pca_F$Sites)
pca_F$Sites=ifelse(pca_F$V1%in%SNI$V1,"SNI",pca_F$Sites)

cols <- c("#0072B2","#009E73","#E69F00")

### Global plotting controls
dataset_name <- "Fast"
jitter_width  <- 0.05
jitter_height <- 0.05
alpha_val     <- 0.6
point_size    <- 4

### Function to plot any two PCs and save to file
plot_pca <- function(df, eigenval, pc_x, pc_y, cols, outdir = ".", format = "pdf") {
  
  # Build the ggplot object
  p <- ggplot(df, aes_string(x = paste0("V", pc_x+2), y = paste0("V", pc_y+2), colour = "Sites")) +
    geom_jitter(width = jitter_width, height = jitter_height, alpha = alpha_val, size = point_size) +
    scale_color_manual(values = cols) +
    geom_hline(yintercept = 0) +
    geom_vline(xintercept = 0) +
    theme_bw() +
    labs(
      x = paste0("PC", pc_x, " (", round(100 * eigenval[pc_x,1] / sum(eigenval[,1]), 2), "%)"),
      y = paste0("PC", pc_y, " (", round(100 * eigenval[pc_y,1] / sum(eigenval[,1]), 2), "%)")
    ) +
    theme(text = element_text(size = 10))
  
  # Construct filename
  filename <- file.path(outdir, paste0(dataset_name,"_PC", pc_x, "_PC", pc_y, ".", format))
  
  # Save depending on chosen format
  if (format == "pdf") {
    ggsave(filename, plot = p, device = "pdf", width = 6, height = 5)
  } else if (format == "png") {
    ggsave(filename, plot = p, device = "png", width = 6, height = 5, dpi = 300)
  } else {
    warning("Unsupported format. Use 'pdf' or 'png'.")
  }
  
  # Also return plot to R session
  return(p)
}

### plots
# Save to PDF (default) in current working directory
plot_pca(pca_F, eigenval, 1, 2, cols)    # PC1 vs PC2

plot_pca(pca_F, eigenval, 1, 3, cols)    # PC1 vs PC3

#plot_pca(pca, eigenval, 1, 4, cols)    # PC1 vs PC4
#plot_pca(pca, eigenval, 1, 5, cols)    # PC1 vs PC5
#plot_pca(pca, eigenval, 1, 6, cols)    # PC1 vs PC6
plot_pca(pca_F, eigenval, 2, 3, cols)    # PC2 vs PC3

plot_pca(pca_F, eigenval, 2, 4, cols)    # PC2 vs PC4

#plot_pca(pca, eigenval, 2, 5, cols)    # PC2 vs PC5
#plot_pca(pca, eigenval, 2, 6, cols)    # PC2 vs PC6
plot_pca(pca_F, eigenval, 3, 4, cols)    # PC3 vs PC4

#plot_pca(pca, eigenval, 3, 5, cols)    # PC3 vs PC5
#plot_pca(pca, eigenval, 3, 6, cols)    # PC3 vs PC6


# Save to PNG in a specific folder
#plot_pca(pca, eigenval, 1, 2, cols, outdir = "/Users/mayroberts/Documents/Projects/LOMINS/Analysis/plink/PCA_locin25psampls/plots", format = "png")

OrigEnz only + SNI PCA plot

## FastEnz only + SNI PCA plot
OrigEnz=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/OrigEnz.txt")
SNI=read.table("/Users/mayroberts/Documents/Projects/FastEnztest/plink_PCA/SNI.txt")
pca_O$Sites <- NA  # create empty column first

pca_O$Sites=ifelse(pca_O$V1%in%OrigEnz$V1,"FastEnz",pca_O$Sites)
pca_O$Sites=ifelse(pca_O$V1%in%SNI$V1,"SNI",pca_O$Sites)

cols <- c("#0072B2","#009E73","#E69F00")

### Global plotting controls
dataset_name <- "Fast"
jitter_width  <- 0.05
jitter_height <- 0.05
alpha_val     <- 0.6
point_size    <- 4

### Function to plot any two PCs and save to file
plot_pca <- function(df, eigenval, pc_x, pc_y, cols, outdir = ".", format = "pdf") {
  
  # Build the ggplot object
  p <- ggplot(df, aes_string(x = paste0("V", pc_x+2), y = paste0("V", pc_y+2), colour = "Sites")) +
    geom_jitter(width = jitter_width, height = jitter_height, alpha = alpha_val, size = point_size) +
    scale_color_manual(values = cols) +
    geom_hline(yintercept = 0) +
    geom_vline(xintercept = 0) +
    theme_bw() +
    labs(
      x = paste0("PC", pc_x, " (", round(100 * eigenval[pc_x,1] / sum(eigenval[,1]), 2), "%)"),
      y = paste0("PC", pc_y, " (", round(100 * eigenval[pc_y,1] / sum(eigenval[,1]), 2), "%)")
    ) +
    theme(text = element_text(size = 10))
  
  # Construct filename
  filename <- file.path(outdir, paste0(dataset_name,"_PC", pc_x, "_PC", pc_y, ".", format))
  
  # Save depending on chosen format
  if (format == "pdf") {
    ggsave(filename, plot = p, device = "pdf", width = 6, height = 5)
  } else if (format == "png") {
    ggsave(filename, plot = p, device = "png", width = 6, height = 5, dpi = 300)
  } else {
    warning("Unsupported format. Use 'pdf' or 'png'.")
  }
  
  # Also return plot to R session
  return(p)
}

### plots
# Save to PDF (default) in current working directory
plot_pca(pca_O, eigenval, 1, 2, cols)    # PC1 vs PC2

plot_pca(pca_O, eigenval, 1, 3, cols)    # PC1 vs PC3

#plot_pca(pca, eigenval, 1, 4, cols)    # PC1 vs PC4
#plot_pca(pca, eigenval, 1, 5, cols)    # PC1 vs PC5
#plot_pca(pca, eigenval, 1, 6, cols)    # PC1 vs PC6
plot_pca(pca_O, eigenval, 2, 3, cols)    # PC2 vs PC3

plot_pca(pca_O, eigenval, 2, 4, cols)    # PC2 vs PC4

#plot_pca(pca, eigenval, 2, 5, cols)    # PC2 vs PC5
#plot_pca(pca, eigenval, 2, 6, cols)    # PC2 vs PC6
plot_pca(pca_O, eigenval, 3, 4, cols)    # PC3 vs PC4

#plot_pca(pca, eigenval, 3, 5, cols)    # PC3 vs PC5
#plot_pca(pca, eigenval, 3, 6, cols)    # PC3 vs PC6


# Save to PNG in a specific folder
#plot_pca(pca, eigenval, 1, 2, cols, outdir = "/Users/mayroberts/Documents/Projects/LOMINS/Analysis/plink/PCA_locin25psampls/plots", format = "png")

Population stats

Calculate Nei’D

stamppNeisD(gl_FastEnzTest, pop = TRUE)

##            [,1]     [,2]     [,3]
## GUorig 0.000000 0.372777 0.028816
## SNI    0.372777 0.000000 0.353807
## GUFast 0.028816 0.353807 0.000000

Calculate Fst

stamppFst(gl_FastEnzTest, nboots = 100, percent = 95, nclusters = 1)

## $Fsts
##           GUorig       SNI GUFast
## GUorig        NA        NA     NA
## SNI    0.6230037        NA     NA
## GUFast 0.0142895 0.5996783     NA
## 
## $Pvalues
##        GUorig SNI GUFast
## GUorig     NA  NA     NA
## SNI      0.00  NA     NA
## GUFast   0.02   0     NA
## 
## $Bootstraps
##   Population1 Population2            1            2            3           4
## 1      GUorig         SNI  0.608020749  0.609649287 0.6097027943 0.610013801
## 2      GUorig      GUFast -0.001724113 -0.001111937 0.0009658338 0.001874286
## 3         SNI      GUFast  0.581142733  0.582386482 0.5826943708 0.588506529
##            5           6           7           8           9          10
## 1 0.61082096 0.612698251 0.613026917 0.613793467 0.614021441 0.614168696
## 2 0.00335393 0.003504329 0.003625191 0.004436341 0.005234761 0.005555987
## 3 0.58905292 0.589635219 0.589814130 0.590006282 0.590242205 0.590792421
##            11         12          13          14          15         16
## 1 0.614580542 0.61487796 0.615239167 0.615292896 0.616101623 0.61660828
## 2 0.005710504 0.00579475 0.006220175 0.006251671 0.006864236 0.00694198
## 3 0.590947054 0.59117345 0.591232918 0.591676135 0.592172909 0.59302375
##            17          18        19          20         21          22
## 1 0.616697685 0.617103542 0.6171496 0.617253914 0.61737842 0.617475599
## 2 0.007091702 0.007258867 0.0075095 0.007823758 0.00820538 0.008613906
## 3 0.593073910 0.593104727 0.5932129 0.594027997 0.59423082 0.594254409
##            23         24          25         26          27          28
## 1 0.617790116 0.61811559 0.618207207 0.61872311 0.618754534 0.619016572
## 2 0.009061875 0.00907548 0.009170421 0.00930843 0.009437779 0.009649831
## 3 0.594861648 0.59498917 0.595345779 0.59535378 0.595399065 0.595457291
##            29          30         31         32         33         34
## 1 0.619087477 0.619180607 0.61918163 0.61927831 0.61930298 0.61947546
## 2 0.009650565 0.009796691 0.01030869 0.01034701 0.01046347 0.01052569
## 3 0.595458509 0.595461510 0.59556184 0.59618615 0.59641782 0.59651156
##           35         36         37         38         39         40        41
## 1 0.61978898 0.61979396 0.62026724 0.62041468 0.62067817 0.62075168 0.6208859
## 2 0.01087547 0.01091272 0.01120417 0.01138633 0.01177514 0.01184398 0.0119559
## 3 0.59670675 0.59689920 0.59701669 0.59718156 0.59750958 0.59754377 0.5975673
##           42         43         44         45         46         47         48
## 1 0.62113415 0.62115287 0.62124665 0.62147058 0.62172748 0.62179902 0.62200498
## 2 0.01209523 0.01219044 0.01224498 0.01230399 0.01232285 0.01244308 0.01246728
## 3 0.59766864 0.59777688 0.59788358 0.59799762 0.59825695 0.59825751 0.59861867
##           49         50         51         52         53        54         55
## 1 0.62243826 0.62279479 0.62283002 0.62291060 0.62299208 0.6232466 0.62346287
## 2 0.01253975 0.01258071 0.01281685 0.01283838 0.01284931 0.0128697 0.01300003
## 3 0.59892088 0.59899943 0.59906284 0.59921410 0.59951879 0.5997534 0.60020506
##           56         57         58         59         60         61         62
## 1 0.62388796 0.62397485 0.62403190 0.62407573 0.62488128 0.62523938 0.62529706
## 2 0.01371166 0.01373416 0.01389332 0.01434193 0.01439945 0.01441518 0.01490872
## 3 0.60076469 0.60116439 0.60135924 0.60161180 0.60173958 0.60186667 0.60238484
##           63        64         65         66         67         68         69
## 1 0.62544380 0.6255104 0.62572758 0.62581673 0.62638348 0.62671957 0.62675402
## 2 0.01514739 0.0152280 0.01526081 0.01539994 0.01612105 0.01637586 0.01676796
## 3 0.60255780 0.6027125 0.60281105 0.60303935 0.60304040 0.60337762 0.60350005
##           70         71         72        73         74         75         76
## 1 0.62681336 0.62703235 0.62708198 0.6271434 0.62746204 0.62757096 0.62768502
## 2 0.01678145 0.01691997 0.01722414 0.0173890 0.01750387 0.01750908 0.01781178
## 3 0.60402763 0.60419690 0.60420446 0.6044011 0.60456610 0.60493709 0.60514837
##           77         78        79         80         81         82         83
## 1 0.62770232 0.62772267 0.6278493 0.62801754 0.62809910 0.62829406 0.62890897
## 2 0.01789464 0.01792016 0.0180338 0.01850626 0.01853681 0.01868865 0.01873624
## 3 0.60523874 0.60525971 0.6053732 0.60542826 0.60556868 0.60559283 0.60565109
##           84        85         86         87         88         89         90
## 1 0.62933035 0.6294674 0.62964218 0.62982907 0.62991453 0.63121758 0.63187351
## 2 0.01917819 0.0192646 0.01933373 0.01941165 0.02010527 0.02024507 0.02042847
## 3 0.60666831 0.6067654 0.60700967 0.60742046 0.60768494 0.60974761 0.61002143
##           91         92         93         94         95         96         97
## 1 0.63205622 0.63243393 0.63299816 0.63444026 0.63457478 0.63459144 0.63589101
## 2 0.02047936 0.02055784 0.02057224 0.02065127 0.02090403 0.02112547 0.02145911
## 3 0.61014876 0.61044159 0.61132850 0.61219420 0.61245124 0.61256387 0.61314668
##           98         99       100 Lower bound CI limit Upper bound CI limit
## 1 0.63731090 0.63744442 0.6386691         0.6097027943           0.63589101
## 2 0.02211402 0.02324004 0.0252369         0.0009658338           0.02145911
## 3 0.61362154 0.61381994 0.6138233         0.5826943708           0.61314668
##   p-value       Fst
## 1    0.00 0.6230037
## 2    0.02 0.0142895
## 3    0.00 0.5996783